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Abstract. The directed-loop scheme is a framework for generalized loop-type updates in quan- 
tum Monte Carlo, applicable both to world-line and stochastic series expansion methods. Here, the 
directed-loop equations, the solution of which gives the probabilities of the various loop-building 
steps, are discussed in the context of the anisotropic 5=1/2 Heisenberg model in a uniform mag- 
netic field. This example shows how the directed-loop concept emerges as a natural generalization 
of the conventional loop algorithm, where the loops are selfavoiding, to cases where selfintersection 
must be allowed in order to satisfy detailed balance. 



1. INTRODUCTION 

Loop algorithms [1, 2, 3] have dramatically improved the performance of world-line 
quantum Monte Carlo calculations [4]. The autocorrelation times can be reduced by 
several orders of magnitude relative to standard local updating schemes [5]. However, 
the conventional loop updates are restricted to certain models and/or limited regions of 
their parameter spaces. In particular, external fields cannot be taken into account when 
constructing a loop, and the loop-flip is then conditional upon a subsequent Metropolis 
[6] accept/reject step. The acceptance probability for large loops in a high field is small, 
and this approach is therefore feasible only at high temperatures or very weak fields [7]. 
The restriction is analogous to that in classical Monte Carlo, where cluster algorithms 
[8, 9] also are not applicable to spin models in a magnetic field. Remarkably, two recent 
generalizations of the loop concept have overcome this problem for quantum systems. 
The worm algorithm [10] for world-lines in continuous imaginary time and the operator- 
loop algorithm [11] for stochastic series expansion (SSE) [12] generalize the loop by 
allowing it to selfintersect and backtrack. The original prerequisite of a cluster algorithm, 
i.e., to express the partition function using new auxiliary variables [8, 13], is then 
circumvented, and the generalized loops can therefore take complicated interactions and 
external fields into account. The loop-building takes place in an extended configuration 
space of the original variables (spin states or occupation numbers), where configurations 
with uncompleted loops (or worms) do not contribute to the partition function (they 
correspond to violation of a conservation law). Such a method was attempted already in 
the early days of the world-line algorithm [14], but without enforcing detailed balance in 
the loop construction. Due to the low acceptance probability for random-walk loops, the 
method was not as efficient as simple local updating schemes [4, 15]. In the worm and 
SSE operator-loop algorithms, detailed balance is ensured by local probabilistic rules, 
and the resulting closed-loop configurations are always accepted. 
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The exact relationship between the conventional loop algorithms [1,2,5] and the more 
general loop-type algorithms allowing selfintersection and backtracking [10, 11] was not 
immediately clear. In particular, the general algorithms did not reduce to a standard loop 
algorithm in regions of parameter space where such an algorithm could be applied. Al- 
though the efficiency was dramatically improved over local updates [16, 11, 17, 18], 
the standard loop algorithm was still much more efficient when applicable. Particularly 
unsatisfying was the fact that simulations could not be carried out as effectively close to 
a region of an applicable conventional loop algorithm as within such a region. This "al- 
gorithmic discontinuity" problem was solved with the introduction of the directed-loop 
algorithm [19], which can often be tuned so that the probabilities of selfintersection and 
backtracking smoothly vanish as a region of an applicable loop algorithm is approached. 
The directed-loop algorithm then becomes identical to a standard single-loop algorithm 
(i.e., one loop at a time is constructed, as in the classical Wolff cluster algorithm [9]). The 
directed loops thus emerge as a natural generalization of the original [1] loop concept, in 
a way similar to the Kandel-Domany generalization [20] of classical cluster algorithms. 

In the directed-loop scheme, the detailed-balance conditions lead to a set of coupled 
equations for the probabilities of the various loop-building (or worm 1 ) steps. These 
directed-loop equations often have an infinite number of solutions, which hence should 
be optimized. The directed-loop algorithm was first developed for SSE, but adaptations 
to world-lines, both in discrete and continuous imaginary time, were also presented in 
the same article [19]. Conceptually, the scheme is simpler (and often more efficient) 
for SSE, and here it will therefore be discussed only within this representation. For 
simplicity, only the S = 1/2 XXZ model (the anisotropic Heisenberg model in a uniform 
magnetic field) will be considered. In this case the optimization criterion for the directed 
loops is taken to be the minimization of the backtracking probability. 

In Sec. 2 the basics of the SSE method are reviewed, first in general and then focusing 
on the details for the S = 1/2 XXZ model. The structure of the operator loops and 
the derivation of the directed-loop equations are discussed in Sec. 3. Some recent 
applications and extensions of the directed-loop algorithm are summarized in Sec. 4. 



2. STOCHASTIC SERIES EXPANSION 

The SSE method [12] is an efficient and widely applicable generalization of Hand- 
scomb's [21] power-series method. To construct the SSE representation of the partition 
function, Z = Tr{exp(— /3//)}, the Hamiltonian is first written as a sum, 

H = ~Y^H a ^ (1) 

a b 

where in a chosen basis (|oc)} the operators satisfy H a ^ b \a) ~ \oc'), where \ a) and |a ; ) 
are both basis states. The subscripts a and b refer to the operator types (various diagonal 



Generally speaking, a worm is a different name for an incomplete loop, but it should be noted that the 
worm-building processes in the worm algorithm [10] differ from those used in SSE operator loops [11, 19] 
and the directed loops for world-lines in continuous or discrete imaginary time [19]. 
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and off-diagonal terms) and the lattice units over which the interactions are summed 
(e.g., the bonds corresponding to two-body interactions). A unit operator Hq q = 1 is also 
defined. Using the Taylor expansion of exp(— j5H) truncated at order M, the partition 
function can then be written as [12] 



a Sm ^ ' 



M 

n // «» 



a), (2) 



where Sm = [a\,bi\, [02,^2]? ••• , [«m>^m] corresponds to the operator product, and n 
denotes the number of non- [0, 0] elements (i.e., the actual expansion-order of the terms). 
M can be adjusted during the equilibration of the simulation, so that it always exceeds 
the highest power n reached; M — > An max , where, e.g., A = 4/3. Then M ~ j5N, where 
N is the number of sites, and the remaining truncation error is completely negligible. 
Defining a normalized state \ ct(p)) as |oc) propagated by the first p operators, 

\a{p))~\[H aiM \a), (3) 

!=1 

the periodicity \a{M)) = |ce(0)) is required for a non-zero contribution to Z. In an SSE 
simulation, transitions (a, Sm) — ► (oc',S' M ) satisfying detailed balance are carried out to 
sample the configurations. Three classes of updates are typically used: 

(i) Diagonal update, where the expansion order n is changed by replacing a fill-in 
unit operator by a diagonal operator from the sum (1), or vice versa, i.e., [0, 0] <-> [d, b], 
where the type-index d corresponds to a diagonal operator in the basis used. 

(ii) Off-diagonal update, where a set of operators {[«p,£>p]} is updated by changing 
only the type-indices a p . Off-diagonal operators cannot be added and removed one-by- 
one with the periodicity constraint \a{M)) = |«(0)) maintained. Local updates involv- 
ing two simultaneously replaced operators can be used [12], but much more efficient 
loop [11] and "quantum-cluster" [22] updates have also been developed. 

(iii) State update, which affects only the state \a) in (2). This state is just one out 
of the whole cycle of propagated states \a(p)), and it can change in the off-diagonal 
updates (ii). However, at high temperatures many sites will frequently have no operators 
acting on them, and they will then not be affected by off-diagonal updates. The states at 
these sites can then instead be randomly modified, as they do not affect the weight. 

Turning now to the anisotropic Heisenberg antiferromagnet in a magnetic field, 

H = jJ^[SjS x J + S y i S y j + ASIS$-hJ^St, (7>0,A>0), (4) 

the standard z-component basis is used: \a) = ...,S Z N ), S^ = ±1/2. Diagonal and 
off-diagonal bond operators are defined, 

Hij, = E + A/4 + h b -AS z i{b) S z m + h b [S z i{b) +S^ (5) 
H 2 , b = -j[Sl {b) S m +S7 {b) S+ b) ], (6) 

where i(b),j(b) are the sites connected by bond b and is the bond-field (e.g., on a d- 
dimensional cubic lattice hi, = h/2d). The Hamiltonian can now be written in the form 
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(1) with a = 1 , 2 and b = 1 , . . .Nb, where Nb is the number of bonds (e.g., Nb = dN on a 
J-dimensional cubic lattice). Note again that the unit operator Hq q = I is not part of the 
Hamiltonian; it is a fill-in element for augmenting the products of order n < M in (2). 

The constant e + ^ + hb has been added to the diagonal bond-operator (5) in order to 
render all its matrix elements positive (e > 0). On a bipartite lattice, the minus-sign in 
the off-diagonal operator (6) is irrelevant, and the expansion (2) is then positive-definite. 
The sign problem for frustrated Xy-interactions [23] will not be considered here. 

Storing the operator sequence Sm and a single state \ oc(p)) (initially \a) = |a(0))), di- 
agonal updates of the form [0, 0] <-> [l,b] can be carried out sequentially for p = 1,...,M 
at all elements [a p , b p ] in Sm with a p = 0, 1. The Metropolis acceptance probabilities for 
such substitutions are [12] 

P([0,0]-[l,fc]) = NbP(S]{p)S){p)\H hb \S]{p)S){p))/{M-n), (7) 
P([l,&]->[0,0]) = (M-n + \)/[N b p{S\{p)S){p)\H l)b \S]{p)S){p))], (8) 

where, as always [6], P > 1 should be interpreted as probability one. The spins S$(p) 
refer to the propagated states (3), which are generated one-by one during the diagonal 
update by flipping spins whenever off-diagonal operators [2,b] are encountered. 

In the early applications of the SSE scheme [12], local off-diagonal updates involving 
simultaneous substitution of two operators were used, i.e., [l,fc p ][l,£> 9 ] <-> [2,b p ][2,b q ]. 
The operator-loop update [11] to be discussed next is a much more efficient way of 
sampling the off-diagonal operators. 



3. OPERATOR LOOPS AND DIRECTED LOOPS 

To begin the discussion of SSE loop-type updates, it is useful to first consider one of the 
simplest cases; the isotropic model, with A = 1 , h = in (4). In this case, setting e = in 
Eq. (5), both the diagonal ([1, b]) and off-diagonal ([2, b]) operators can act only on anti- 
parallel spins, and the corresponding matrix elements are 1/2 [neglecting the negative 
sign in (6)]. Fig. 1 shows a graphical representation of a valid configuration, along with 
an illustration of a deterministic loop update [11]. Here all fill-in operators [0,0] and the 
corresponding propagated states have been left out since they are irrelevant in the loop 
update (as the expansion-order n does not change). Selecting the starting point and a 
direction (up or down) at random, a loop is constructed using a completely deterministic 
rule: Moving along the chosen direction, whenever an operator is encountered the path 
switches to the other spin connected to that operator, and the direction of movement 
is reversed. This will eventually lead to a closed loop when the initial starting point 
is reached. A new valid configuration is then obtained by flipping the spins along the 
loop and changing the types of all operators encountered; diagonal <-> off-diagonal (in 
practice, the changes are carried out on the run while building the loop). Operators 
encountered twice will remain unchanged. Since all non-zero matrix elements of the 
bond operators equal 1/2, the new configuration has exactly the same weight as the 
old one, and the loop-flip can hence always be accepted. This type of loop is self- 
avoiding by construction. Instead of constructing loops one-by-one and flipping them 
with probability one, the configuration can therefore also be decomposed into all its 
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FIGURE 1. An SSE configuration of order n = 10 for a 6-site isotropic XXZ chain. Open and solid 
circles correspond to up and down spins, and open and solid bars represent diagonal and off-diagonal 
operators, respectively. The construction of a loop is illustrated to the left, where the start/end point is 
indicated with a bar/arrow, and he initial direction of movement is upward. The configuration obtained 
when the loop has been flipped is shown to the right. 



loops (each spin belongs to exactly one loop), which are flipped independently of each 
other with probability 1/2 (as in the classical Swendesen-Wang [8] algorithm). 

The deterministic loop update clearly relies on the symmetry of the model at the 
isotropic point (A = 1, h = 0). The simple rule of "switch and reverse" when an operator 
is encountered has to be modified in order to construct a more general loop-scheme, ap- 
plicable for any A, h. In the general operator-loop update [11], there are four possibilities 
for the worm-like path to proceed when an operator is encountered; it can continue on 
the same spin or switch to the other spin connected by the operator, and in either case 
the direction of the movement can be up or down. The directed-loop approach provides 
the general detailed-balance conditions that these probabilities have to satisfy. 

In order to discuss the general operator-loop update and the directed-loop scheme, 
it is useful to introduce a different representation of the SSE configurations. It is not 
necessary to store the full states \ct(p)) shown in Fig. 1; the same-spin "lines" between 
the operators clearly contain a great deal of redundant information. One can represent 
the matrix element in Eq. (2) as a linked lists of vertices [11]. Note first that the weight 
of a configuration (a, 5m) can be written as 

W(a,S M ) = P n{M m n)l flW(p), (9) 

where the product is over the n non-[0,0] operators in Sm- W(p) will be referred to as 
a bare vertex weight; it can be written as a matrix element of the full bond operator 
H b = Hi b + H 2j b at position p; 

w(p) = (^)(p)^)(p)I%I%)(p- ^ z j{bp) (p- 1)>- do) 

A vertex represents the spins on bond b p before and after the operator has acted. These 
four spins constitute the legs of the vertex. There are six allowed vertices, with four 
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FIGURE 2. (a) All vertices for the S = 1 /2 XXZ model, with their vertex weights W k . (b) The linked- 
vertex representation (right) of a full 3-spin SSE configuration with n = 3 (left). 



different vertex-weights as illustrated in Fig. 2(a). The weights are 

wi = mm id = e , 

W 2 = (IT \H b \ IT) = W 3 = (TII^ITI) = A/2 + ^ + e, (11) 

w 4 = (TII^IIT) = w 5 = (ITI^ITI) = 1/2, 

w 6 = (TT \H b \ TT) = 

An example of a linked-vertex representation of a term with three bond operators is 
shown in Fig. 2(b). The links connect vertex-legs on the same site, so that from each leg 
of each vertex, one can reach the next or previous vertex-leg on the same site (i.e., the 
links are bidirectional). In cases where there is only one operator acting on a given site, 
the corresponding "before" and "after" legs of the same vertex are linked to each other 
[as is the case with the legs on site 1 in Fig. 2(b)]. 

The building of a loop in the linked-vertex representation consists of a series of steps, 
in each of which a vertex is entered at one leg (the entrance leg) and an exit leg is chosen 
according to probabilities that depend on the entrance leg and the spin states at all the 
legs [i.e., the vertex type, k = 1, . . . , 6, in Fig. 2(a)]. The entrance to the following vertex 
is given by the link from the chosen exit leg. The spins at all visited legs are flipped, 
except in the case of a bounce, where the exit is the same as the entrance leg, and 
only the direction of movement is reversed. The starting point of the loop is chosen 
at random. Two link-discontinuities (which are analogous to the source operators in 
the worm algorithm [10]) are then created when the first entrance and exit spins are 
flipped, i.e., these legs will now be linked to legs with different spins . Configurations 
contributing to Z only contain links between same-spin legs [as in Fig. 2(b)]. When 
the loop closes, the two discontinuities annihilate each other, and a new contributing 
configuration has then been generated. 

The probabilities for the different exit legs (e = 1, . . .4), given the type of the vertex 
(k = 1,...6) and an entrance leg (i = 1,...4), are chosen such that detailed balance 
is satisfied. This leads to the directed-loop equations, which are constructed in the 
following way: Unknown weights a e (i : k) are first assigned to all possible paths (z — > e) 
through each vertex k. The sum of all these path weights over all exits e must equal 
the bare vertex weight (10), i.e., the matrix element before the entrance and exits spins 
have been flipped; Y* e a e(hk) = W^. The actual normalized exit probability is the path 
weight divided by the bare vertex weight; P e (i,k) = a e (i : k) /W^. The key observation 
leading to the directed-loop equations [19] is that the weights for vertex-paths i — > e that 
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FIGURE 3. Two closed sets of vertex paths, with their corresponding bare vertex weights Wt and path 
weights dj, bj. The entrance legs are indicated with arrows pointing into the bare vertices, and the exit 
legs for the three allowed paths are at the arrows pointing out from the vertices. The entrance and exit 
spins on the paths have been flipped. 



constitute each other's reverses have to be equal: If the path i — > e through vertex k leads 
to the vertex k' when the entrance and exit spins have been flipped, then the reverse path 
e — > / through k' yields vertex k, and if a e (i,k) = a,i{e,k') it immediately follows that 
P e (i, k)Wt = Pi(e 7 k')Wtf, i.e., local detailed balance is satisfied. 

The condition a e {i,k) = ai{e,k') couples some of the equations Y,e a e(i,k) — f° r 
different entrance legs i and vertices k, and these equations have to be solved for all 
the weights a e {i,k). Typically, not all the equations are coupled, however, but there are 
several different sets that can be solved independently of each other. Such closed sets for 
the XXZ model are illustrated in Fig. 3. Here the path- weights are labeled <2 ; and bi for 
the two sets (a),(b), and paths that constitute each other's reverses have been assigned 
the same weight. Note that one of the four exits always leads to a new vertex that does 
not correspond to a term in the XXZ Hamiltonian; these paths are not allowed and are 
not included in the figure. The directed-loop equations for the two closed sets are 

W\ = a\ +a 2 + a 4 , W 6 = b\+b2 + b4, 

W 3 = a 1 +a 3 +a 5 , W 2 = b 1 +b 3 + b 5 , (12) 

W 4 = a 2 + a 3 +a 6 , W 5 = b 2 + b 3 + b 6 , 

where the bare vertex weights are given in Eq. (11). All the remaining closed sets 
are related to those in Fig. 3 by trivial symmetries, and for h = the corresponding two 
sets of equations (12) are identical. Note that there are six weights ai and bi to be solved 
for in each set, but only three equations. There is thus an infinite number of solutions, 
even with the requirement that all weights have to be positive (since the probabilities are 
obtained by dividing by the positive matrix elements W^). 

In Ref. [11], a particular "heat-bath" solution was obtained by working directly with 
the probabilities, instead of analyzing the path-weights of the directed-loop scheme. The 
directed-loop equations (12) provide a more general framework for finding the optimal 
solution, i.e., the one which leads to simulations with the shortest autocorrelation times. 
There is currently no rigorous way of finding the optimal solution, but heuristic argu- 
ments have been put forward [19]: It is a reasonable assumption that the probabilities for 
the bounce processes (i.e., the last columns in the sets in Fig. 3) should be minimized, 
as they do not accomplish any vertex changes and cause the loop-building process to 
backtrack one step (and sometimes more than one step as the loop-building continues 
in the opposite direction). For the model considered here, minimizing the bounces leads 
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FIGURE 4. (a) The shaded areas show the regions of the A > Q,h > parameter space where the 
solutions of the directed-loop equations are bounce-free (set a left, set b right), (b) Algorithmic phase 
diagram. The labels A,B,C,D correspond to the solutions for the weights listed in Table. 1 . 



to a unique solution. Both sets (a) and (b) have regions where all bounce probabilities 
are zero, as shown in Fig. 4(a). The "algorithmic phase" diagram for the full minimum- 
bounce solution has four different regions, as shown in Fig. 4(b). They correspond to 
different analytical forms of the solution. All the path-weights in these regions are listed 
in Table 1 . As discussed above, the actual probabilities are simply obtained by normal- 
izing with the matrix elements of the corresponding equations in (12). 

The solutions in the regions A-D in Fig. 4(b) are continuous across the boundaries. 
In particular, it can be noted that when the isotropic Heisenberg point, A = 1, hf, = 0, is 
approached, and the minimum value £ m i n of the constant £ is used, the only surviving 
vertex process is the switch-and-reverse, corresponding to the weights a^, in Fig. 3. 
This is exactly the process used in the deterministic update illustrated in Fig. 1. Hence, 
the general directed loops indeed smoothly reduce to these very special symmetric loops, 
that were constructed in a different manner by using the fact that the diagonal and off- 
diagonal matrix elements have the same values 0, 1/2 at the isotropic point. 

It is interesting to note that the constant £ appears only in the weights a\, b\, corre- 
sponding to the "continue- straight" process in Fig. 3. One can always choose £ = e m [ n 
(which has the advantage that it minimizes the average expansion-order (n) and the cut- 
off M), but in some cases when £ m j n = a non-zero £ can lead to shorter autocorrelation 
times [19]. In a world-line formulation of the directed-loop update [19], the constant £ 
does not appear at all, but apart from this the path- weights are the same as in Table 1. 
When taking the continuum limit of the time-discretized world-lines, it has been shown 
[19] that the vertex-probabilities reduce exactly to those obtained before [3, 24] on the 
line h = 0, A < 1 . This smooth reduction of the directed loops to the conventional self- 
avoiding loops shows that this scheme is a natural generalization of the loop concept 
when the requirement of self-avoidance is relaxed. In the earlier generalized loop al- 
gorithms, i.e., the worm algorithm [10] and the SSE operator-loop algorithm with the 
simple heat-bath probabilities [11] (which also correspond to a solution of the directed- 
loop equations), the bounce probability does not vanish as h — > 0, and hence any formal 
relationship between the generalized and conventional loop algorithms was unclear. The 
relationships between the various loop-type methods are also reviewed elsewhere in this 
Volume [25]. 
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TABLE 1. Vertex-path weights and the minimum value of the constant e in the different 
regions of Fig. 4(b). A short-hand notation A ± = (1 ± A)/4, / = hb/2 is used. 
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Practical implementation details of the directed-loop algorithm have not been dis- 
cussed here; they have been outlined in Ref. [19]. Some example-programs are also 
available on-line [26] . 



4. APPLICATIONS AND GENERALIZATIONS 

The directed-loop algorithm has already been applied to several models in addition 
to the 5=1/2 antiferromagnetic XXZ model discussed here and in Ref. [19] (the 
ferromagnetic case can be treated in a very similar manner [19]). An identical algorithm 
in the continuous-time world-line representation has been applied in a large-scale study 
of the weakly anisotropic (A > 1,^ = 0) system [27]. Extensions to higher spins and 
softcore boson models have been explored by several groups [28, 29, 30, 31]. Four-spin 
interactions have also been considered [32]. An application to the ID extended Hubbard 
model has produced high-precision results [33] for larger systems sizes than what was 
practically feasible with previous methods. A different type of directed loops have been 
developed for quantum-rotor models [34]. 

In general, for S > 1/2 and softcore bosons, minimization of the bounce probability 
does not lead to a unique solution of the directed-loop equations [29], and therefore 
some other constraints have to be applied as well. It has also been pointed out that 
the minimization of the bounce is not necessarily the optimal strategy [31] (which was 
also anticipated not to be strictly true from the outset [19]). Also, for the Heisenberg 
model with S > 1, in order to eliminate the bounces completely one has to assign 
multiplicative weights ^ 1 also to the discontinuities (sources) that exist while the loop 
is being constructed [31]. 

It appears that in most cases it is relatively easy to find low-bounce solutions that work 
well in practice, but it remains a challenging problem to find a scheme for automatically 
generating the optimal vertex-probabilities, i.e., without carrying out time-consuming 
tests of actual simulation programs. The strength of the directed-loop scheme is that it 
provides a well-defined, general mathematical framework for this pursuit. 
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